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Abstract: A reduced-order model based on Proper Orthogonal Decomposition (POD) is pro- 
posed for the bidomain equations of cardiac electrophysiology. Its accuracy is assessed through 
electrocardiograms in various configurations, including myocardium infarctions and long-time sim- 
ulations. We show in particular that a restitution curve can efficiently be approximated by this 
approach. The reduced-order model is then used in an inverse problem solved by an evolutionary 
algorithm. Some attempts are presented to identify ionic parameters and infarction locations from 
synthetic ECGs. 

Key-words: Cardiac electrophysiology, reduced-order model, inverse problem, POD 
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Modeles reduits pour I'electrophysiologie 
cardiaque. Application a I'identification des 

parametres 

Resume : Un modele reduit base sur la decomposition orthogonale aux 

valours proprcs (POD) est propose pour les equations bidoniaines de I'electro- 
physiologie cardiaque. Sa precision est evaluee a travers des electrocardio- 
grammes dans diverses configurations (modeles d'infarctus du myocarde, si- 
mulations de longue durce, etc.). On montre en particulier que la courbe de 
restitution pent etre approchee de maniere efficace par cette approche. Le mo- 
dele reduit est ensuite utilise dans un probleme inverse resolu par un algorithme 
evolutionnaire. On prcsente des rcsultats d 'identification de parametres ioniques 
et de localisation d'une zone infarcie a partir d'un EGG synthetique. 

Mots-cles : Electrophysiologic cardiaque, modeles reduits, probeme inverse, 
POD 
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1 Introduction 

The bidomain equations used to model the cardiac electrophysiology is known 
to be very demanding from a computation viewpoint j8j [29] . To reduce its com- 
plexity, reduced-order models based on physical arguments can be considered. 
For example, assuming that the anisotropy ratios of the electrical conductivity 
tensors are the same in the intra and extracellular media, a simplified model 
called "monodomain" can be derived [3 IH] . Another example of physical simpli- 
fication is given by the Eikonal model that essentially describes the propagation 
of the activation front [TBI Here we follow another route: keeping the 

physics of the model, we reduce the complexity of its discretization by using a 
Galerkin basis adapted to the kind of solutions that are looked for. This basis 
can be computed by different approaches. In this work, we choose the most 
popular one, which is known as the Proper Orthogonal Decomposition (POD) 
or Karhunen-Loeve method. This approach has been used in many fields of 
science and engineering, but to the best of our knowledge not for the equations 
of cardiac electrophysiology. The present study is a first step in this direction. 
We consider various configurations of interest in order to point out those where 
the reduced-order model seems promising and those where it has to be used 
with much care. 

After a brief presentation of the equations and the numerical methods used 
for the full and reduced-order methods in Section |2] the reduced-order model 
is tested in three situations in Section [3] perturbation of some parameters, 
long-time simulation and infarction modeling. In all these cases, the results are 
assessed through the corresponding electrocardiogram (EGG). In Section |4] the 
reduced-order model is used with an evolutionary algorithm to identify some 
parameters of the problem and the location of infarcted regions from the elec- 
trocardiograms. The main conclusions of the study are then summarized in 
Section [5l 



2 Models and methods 
2.1 Models 

The electrical activity in the heart is modeled by the bidomain equations (see 
pn [251 [551 [3lJ e.g.). The extracellular potential Ue and the transmembrane 
potential Vm are solutions in the heart domain flu of the following system 

(C'm^ + /io„(V„i, W)^ - div(criVK„) - div((TiVUe) = ^m4pp m (0, T) X 

- div(((Ti + o-e) Vuc) - div(criVKi) = in (0, T) x 

(1) 

where is the rate of membrane area per volume unit, the membrane 
capacitance per area unit, cr; and cTc the intra- and extracellular conductivity 
tensors. 

The term /app is a given source function, used in particular to initiate the 
activation, and lioniVrmw) represents the ionic current across the membrane. 
The ionic variable w is solution of the ordinary differential equation 

^4-5(K.,u.) = 0in(0,r)xr!H. (2) 
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In this study, the dynamics of w and /jon are described by the phenonieno- 
logical two- variable model introduced by Mitchell and SchaefFer [20]: 



W 1 

2 if < V^ate, 



i(^'^ax ^nin) 
W 



(3) 



'^closc 



if > Vgate, 



where Tin, Tout, Topen, Tdosc, ^gato are given parameters and Kiin, Vmax are scal- 
ing constants (typically -80 and 20 mV, respectively). Although this model de- 
scribes in a very simplified way the ionic exchanges, it allows to recover realistic 
action potentials at the cell scale. Moreover, each parameter has a physiologi- 
cal meaning: the time constants ri„. Tout are respectively related to the length 
of the depolarization and repolarization (final stage) phases, Topcn and Tdosc 
are the characteristic times of gate opening and closing respectively and V^ato 
corresponds to the change-over voltage. 

On the boundary of the heart domain, it is generally admitted that the 
intracellular medium is isolated. So we have 

(Ti'VV^aL ■ n + a-iVuc ■ n = on SI^h- (4) 

For the second boundary condition, two kinds of coupling between the heart and 
the external medium - that will be called the torso - are usually considered. The 
most complex one is to assume that the potentials and currents are continuous 
between the extracellular and the torso, which corresponds to a strong coupling. 
In the present study, the effect of the torso on the heart is neglected and it is 
assumed that the extracellular medium is isolated (see [SKS]). This corresponds 
to a weak coupling described by the following conditions: 

on dflii, 
on oilii. 

In the torso domain, denoted by ilx, the electrical potential ut is solution 
of the generalized Laplace equation: 

div(crxVuT) — 0, in f2x- (6) 

where ctt stands for the conductivity of the different tissues (lungs, bones, etc.). 
Assuming that the body surface Text is isolated, we have 

cttVut • n = 0, on Text- (7) 

This model allows to reproduce the main recording of the electrical cardiac 
activity made by medical doctor, the electrocardiogram (EGG). It is a set of 
12-graphs of voltages vs. time obtained from electrodes located on standard 
positions of the surface of the body. 

In [3J , under appropriate modeling assumptions (anisotropy, cell heterogene- 
ity, simplified model of His bundle) , realistic EGGs are provided in the healthy 
case and for two pathologies: left bundle branch block and right bundle branch 
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block. The ECGs obtained satisfy the usual criteria used by medical doctors 
to detect these pathologies. We refer to [3_ for a description of the modeling 
choices made in our study. 

The heart-torso uncoupling conditions (|5| induces a significant reduction of 
the computational cost for the numerical resolution since the heart problem can 
be solved independently of the torso problem (in our case, the weak coupling 
is about 60 times faster than the strong coupling). As noticed in [3_, making 
this hypothesis increases the amplitude of the EGG but the shape of the EGG is 
only slightly modified. Thus, this simplification will be used in the present study. 



2.2 Numerical methods for the full-order model 

To solve numerically the model given by ([l|-([7|, we follow the procedure ex- 
plained in [2 with the heart-torso uncoupling assumption. The spatial dis- 
cretization is based on the finite element method and the time discretization is 
performed using the second order BDF implicit scheme with an explicit treat- 
ment of the ionic current. 



Let us first introduce some notations to present the resolution algorithm. 
Let Vji C 7?^ (fin) and Wh C H^{ilT) be two finite dimensional subspaces of 
continuous piecewise affine functions. We define Wh,Q — {ipT ^ Wh/ifir — 
on diln}- 

We assume that the time interval [0,r] is divided in K subintervals [ife, ife+i], 

< k < K ~ 1, with tk kSt where 6t T/K is the time step. We 
denote by {V^,Ug,w'^) the approximated solution obtained at time tk- Then, 
(y^"*"^, u;*^"*"-^) is computed as follows: 

• First, a second order extrapolation of T^n at time tk+i and G Vh are 

computed: 

Vj^+^ = 2Vj^ - V^r^ and 



1 

St \ 2 



Then, {V^l+^,u^+^) e V,, x Vh, with / u^'+i = is computed as the 
solution of: 



(8) 



for all (0, Vo) eVhX Vh, with / = 0. 



For other strategies for the time discretization, we refer to [2 |T2 [THl [32 for 
semi-implicit schemes or |15 | l33 l l35] for operator splitting approaches. 
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Due to the second order approximation Vj^'^^ of V^'^^, the proposed time 
scheme has the convenient property to be associated to a constant matrix. In- 
deed, since the ionic current is expHcit, dsl is equivalent to a linear system whose 



matrix in 



j,2N,2N ;„ 



is given by 



(9) 



where matrices M e M^^^, e M^-^ and K2 E M^'^ are given by 

where {4>i)i<i<N is a finite element basis of V/,.. 

Once is known, the torso potential u^^^ e W^;; is obtained by solving: 



(10) 



4+^ = on ai^H- 



Since this last step depends linearly on u^^^ and is weakly coupled to the rest of 
the problem, it can be solved very efficiently by precomputing a transfer matrix 
|3]. Thus, the main computational cost comes from the resolution of system ([s]) 
satisfied by (y^+\u^+i). 



2.3 Reduced-order Modeling 

For completeness, we briefly recall the principle of the the Proper Orthogonal 
Decomposition (POD) method. We refer the reader interested by more details 
to [T71 [23 for example. POD is a method used to derive low-order models 
by projecting the system onto subspaces spanned by a basis of elements that 
contains the main features of the expected solution. To generate the POD basis 
associated with a precomputed solution u ~ {V^,Uc) of the Galerkin problem 
([8| , we make a first numerical simulation (or a set of simulations) and keep some 
snapshots u{tk), I < k < p. Then a singular value decomposition (SVD) of the 
matrix B = {u{ti), . . . ,u{tp)) € E^^'^* is performed: 

B = USV', 

where U G M2JV,2Ar y g j^p.p orthogonal matrices, S £ M^w.p jg ^^e 
matrix of the singular values ordered by decreasing order. The -/Vmodcs first 
POD basis functions {4'i}i<i<A'modos ^^e then given by the iVmodes first columns 
of U and the POD Galerkin problem is solved by looking for a solution of the 
type 

The 2N X 2N sparse matrix given by (|9| is thus replaced by a full matrix of size 
-^modes X -^modcs with the POD method. Since this matrix is constant in time. 
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it is therefore projected on the POD basis and factorized (or even inversed) 
only once at the beginning of the computation. For the simulations presented 
in this paper, the order of magnitude of A'^odcs is typically one hundred and 
the reduced-order model resolution is about one order of magnitude faster than 
the full-order one. 

Note that the model is reduced by changing the basis the solution is approxi- 
mated on, but it still corresponds to the discretization of the original system 
In particular, it depends on exactly the same parameters as the full-order model. 
In general, the reduced-order model does very well at reproducing the solutions 
that have been used to generated the POD basis. But it is difficult to anticipate 
how it behaves when the parameters are modified. This issue is difficult and 
is still the topic of active researches, see for example [1], [6], [14]. In the next 
sections, we will show situations where the reduced model usually works cor- 
rectly and others where it does not. For the latter, we propose some strategies 
to improve the results. 

3 Application to forward problems 

In this section, the reduced-order model is used in different configurations and 
compared to the full-order one. It is important to note that the accuracy is not 
assessed on the whole solution, but only on the EGG, which is considered as the 
"output of interest" of these simulations. 

3.1 POD when some model parameters vary 

A POD basis is computed for a given set of parameters. We consider in this 
section what happens when this basis is used to solve a reduced-order model 
with different parameters. Our investigation is restricted to Tin,Tciosc, and 
Cni because it has been shown in [5] that they are the ones the EGG is the 
most sensitive to. All the simulations are performed on an idealized geometry 
of ventricles based on two ellipsoids (see [30]). The mesh contains 418465 tetra- 
hedra, the time-step is fixed to 0.5 milliseconds and the parameters of the model 
are given in Table [l] 



Am(cm-i) 


Cm(mF) 


Tin 


Tout 


"^opcn 


close 


_cndo 
close 


_mcell 
close 


cpi 
close 


V^ato 


^nin 


^nax 


200 


10--^ 


16 


360 


100 


120 


130 


140 


90 


-67 


-80 


20 



Table 1: Gell membrane parameters. 



Let us first consider a perturbation of Tdosc; which mainly affects the re- 
polarization phase i.e. the T-wave in the EGG. In our model, the heart is 
divided in four regions where Tdoso takes different constant values (see [3]). 
We consider the value in the epicardium of the left ventricle rfj^gp and in the 
right ventricle t^^^,. A POD basis of 80 vectors is first constructed with 
("^ciosc "^ciosc) = (100,100). Then the corresponding reduced-order model is 
solved with (rXc^Tc^oso) = (80,130). Figure [l] (left) shows that the full and 
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reduced models are in good agreement. We only refer to the first lead of the 
ECG for convenience, but the same trend is observed on the other leads. It is 
interesting to note that the experiment used to generate the POD basis (with 
('^cioic'^c^D = (100' 100)) has a negative T-wave (Figure [l] right). It is there- 
fore not trivial that the reduced model is able to provide the correct (positive) 
T-wave. If, instead of considering the ECG only, we compare the extracellular 
potential of the full and reduced models, the relative difference is 2.3 x 10^^, 
in Euclidean norm in Qn x [0,T], with T = 500ms. The full 3D solutions are 
therefore in good agreement too. 




Figure 1: Left: First lead of the ECGs with (rXsc'^dD = (80, 130). Compari- 
son between the full model (dotted red line) and the reduced model (black line) 
with the POD basis generated with (r^Xlc '^dD = (100, 100). Right: First lead 
of the ECG computed with the full-order model for (r^so, = (100, 100). 




time (ms) time (ms) 



Figure 2: Left: First lead of the ECGs with (rin, C^, A^^, r^^^^J = 
(0.8, 10"'^, 200, 120). Comparison between the full model (dotted red line) 
and the reduced model (black line) with a POD basis generated with 
(Tin, Cm, An, rR^,e) = (1-5,2 X 10-3,100,50). Right : First lead of the 
ECG computed with the full-order model for (rin, Cm, ^m, t^osc) = (1-5,2 x 
10-3,100,50) 

Consider now perturbations of Tin, C,n, and t^i^^^- Simulations are run 
with Tin = 0.8, Cm = 10-3, = 200, rX,„ = 120. We compare the ECGs 
obtained with the full-order model and with the reduced-order model corre- 
sponding to a POD basis generated with Tin = 1.5, Cm = 2 x IQ-^, A^ = 
100, T^ia^^ — 50. In spite of some differences (slight temporal shift or differ- 
ence of amplitudes, see Figure |2] left), the results reasonably match. Again, 
this matching is not trivial since the ECG corresponding to the values used to 
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generated the POD basis looks totally different (see Figure |2] right). 

Although it is not possible to mafce a general statement, these results are rep- 
resentative of several experiments that have been performed and not reported 
here: the reduced-order model usually gives an acceptable EGG when the pa- 
rameters TinjTciosc, and Cm vary within a reasonable range. When these 
parameters vary too much, it is recommended to generate new POD bases, clos- 
est to the region of interest. This will be done in Section^ to address inverse 
problems. 

3.2 POD for long-time simulations and restitution curves 

The restitution curve represents the dependence of the action potential dura- 
tion (APD) on the preceding diastolic interval (DI). This curve is known to be 
important in the understanding of some arrhythmias |26l I20| . In addition, it 
has been shown in [19] that the restitution curve can be used to identify some 
parameters in (|3|. The simulation of the restitution curve is extremely chal- 
lenging for 3D models since it requires several dozen of heart beats, whereas the 
simulation of one heart beat is already very demanding. We propose to inves- 
tigate reduced-order modeling to simulate a series of heart beat with variable 
durations. 

We consider a sequence of 12 heart beats, with a heart rate going from about 
55 bpm to about 110 bpm (at each beat, the period is decreased of 50 ms). The 
simulation is first run with the full-order model for 400 milliseconds (less than 
one heart beat) with a timestep of 0.5 milliseconds. A snapshot is taken every 4 
time steps. Then a POD basis of 100 modes is constructed from these snapshots 
and this basis is used for the long-time simulation (10 seconds). For comparison 
purpose, the full-order model was also run for the long-time simulation. Figure|3] 
shows the first lead of the EGG obtained with the full and reduced-order models: 
the two solutions are almost superimposed during the whole simulation. To take 
a closer look, we plot in Figure |4] the corresponding restitution curves obtained 
by recording the APD and DI in a cell localized at the epicardium of the right 
ventricle. The two curves are in very good agreement. The reduced-order model 
therefore seems to be able to correctly handle a heart rate acceleration, even 
when the POD basis has been generated from only one heart beat. Of course, if 
the shape of the propagation front was significantly altered because of the heart 
beat increase, the results would probably deteriorate. But our results show that 
within a significant range of heart rates the reduced-order model is doing very 
well. 

3.3 Full and reduced-order simulation of an infarction 

In the EGG model has been used to simulate healthy cases and bundle 
branch blocks. We propose in the present work to see the effect on the EGG 
of myocardial infarctions. An infarcted region cannot be activated anymore be- 
cause it has been damaged by an interruption of blood supply. In our simplified 
framework, this defect is modeled by dividing by 10 the parameter Tout in the 
infarcted region. Doing so, the region keeps unactivated even when it is stimu- 
lated by the action potential. With more sophisticated ionic models, we could 
for example modify the behavior of the extracellular potassium [25| . 
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First lead : I - first beat First lead : I - last beats 




1000 9000 10000 



Figure 3: First lead of the ECG over a long period (10 seconds) for an increasing 
heart rate (55 beats per minute to 110 beats per minute). On the left-hand side, 
a zoom of the first second of the simulation, on the right-hand side a zoom of the 
last second of the simulation. Solid line indicates the complete model solution, 
dotted line the reduced model one (they are almost superimposed). 



According to electrocardiology books |10| |24] , the main consequence of an 
infarction on a real ECG is an elevation or a depression of the ST segment in dif- 
ferent leads. The magnitude of this elevation or depression and the leads where 
theses changes are visible depend on the position of infarction. Particularly, the 
main features we should find are: 

• in the case of a posterior infarction: a depression in the ST segment in 
the Vi and V2 leads; 

• in the other cases: an elevation in the ST segment with an inverted T 
wave; 

• in the case of an anterior infarction we should look at Vi, V2 or V3, in the 
case of a lateral infarction at 1 or aVL and for an inferior infarction at 
II, III or aVF. 

FiguresJE] and [6] show the simulated ECG for the main kinds of infarction: 
in Figure [Sjthe anterior and the posterior ones, and in Figure [6] the lateral and 
the inferior ones. The result are not very good for the inferior infarction: as 
expected we see an ST elevation in II but we also observe an important ST 
elevation in the last three leads, a depression in III and no sign on aVF which 
is not expected. The difficulty of simulating the inferior infarction is probably 
due to the fact that this zone is very close to the initial activation region. It also 
seems that the QRS complex are not exactly as they should be. Nevertheless, 



Inria 



Reduced-order modeling in cardiac electrophysiology 



11 



Restitution Curve APD AD\ ) 

n+1^ n' 




Dl (msec) 



Figure 4: Restitution curve obtained solving a 12 second simulation with the 
complete model (blue line) and the POD (green line). 




we find as expected a depression (resp. elevation) of the ST segment in Vi, V2 
and V3 leads in presence of a posterior (resp. anterior) infarction. We also find 
a ST elevation in I and aVL in presence of a lateral infarction. In conclusion, 
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the results obtained with the full-order model for the ST segments are very 
satisfactory for the anterior, posterior and lateral infarctions. 

We now propose to investigate the same situations with the reduced-order 
model. Contrary to the conclusions of the previous two sections, we will see 
that a too naive use of POD gives in this case very poor results. 

We first generate the POD basis from a healthy simulation : we run a 400 
milliseconds simulation for an healthy test case, with a time-step of 0.5 mil- 
liseconds, keep snapshots every 2 milliseconds and obtain a basis of 100 vectors. 
Then, we use this basis to simulate an infarction centered in the arbitrary point 
P indicated by the arrow in Figure [9] 

A comparison of the green line (reduced-order model) and blue line (full- 
order model) in Figure |7] shows that this basis does not allow to approximate 
accurately the EGG. 

A look at the transmembrane potential (Figure |8| shows that the "healthy" 
POD basis is indeed unable to capture the sharp variation induced by the in- 
farcted region. This explains the poor results observed on the EGG. 

To improve the space approximation property of the POD basis, we propose 
to collect the snapshots for different infarction points. More precisely, 100 snap- 
shots are taken from an healthy case and 50 snapshots for each infarction of the 
18 points shown in Figure [9] Moreover, since most of the solutions variation 
occurs during the first 100 ms, half of the snapshots are taken in this period and 
the other half between 100 ms and 400 ms. Then a POD basis of 100 vectors is 
computed as usual. 

This new POD basis is used to simulate the infarction on the point P in- 
dicated in Figure |9] A comparison of the results obtained with the full and 
reduced-order models is given in Figure [TT] Although the results look better 
than with a basis coming from a pure healthy case, we observe that the solution 
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Figure 7: Simulated ECGs for an infarction centered in an arbitrary point P 
with a POD basis generated from the healthy case: EGG with the full model 
(blue line), EGG with the POD (green line) and "healthy" full-order EGG (black 
line) . 

Time =100.0 Time =100.0 

(a) Complete model (b) Reduced model 

Figure 8: Simulation of anterior infarction using the complete 
model (left column) and the reduced model (right column) with 
a POD basis generated from the healthy case. 
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seems to superimpose the solutions coming from all the nearest infarction points. 
As a result, the infarcted area seems too large with the reduced-order model 
during the repolarization phase, which induces a difference of the ST elevations 
in the corresponding ECG. Indeed, we see in Figure [TO] that the curves of the 
ECGs have a similar shape but different magnitude for the full and reduced- 
order models. This discrepancy would probably be reduced by refining the grid 
of the precomputed infarcted regions. 



First lead : I 

4 I .. I .... I ..... . 




Figure 9: The mesh and the 18 
points used to build the infarction 
POD basis. The point out of the 
mesh lines, indicated by the pink 
arrow, is the point P considered 
in the example. 



Figure 10: First lead of the ECG for an in- 
farction centered in point P: ECG with the 
complete model (green line) and ECG ob- 
tained using POD (blue line). The healthy 
reference case is given by the black line. 



rime = 70.0 Time = 200.0 



Vl. 

(a) t = 70ms 



v:. 

(b) t = 200ms 



Figure 11: Simulation of an infarction in point P: using the complete model 
(left) and using the POD method (right). 



4 Application to inverse problems 

In this section, we propose a strategy to identify some parameters of the model 
from the ECG. The first application deals with the identification of the 4 pa- 
rameters considered in Section [3.1| In the second application, we try to recover 
the location of an infarction modeled as in Section |3.3[ Note that in this pre- 
liminary study, the reference ECGs used for the identification are "synthetic" 
i.e. previously computed by the model itself. 
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4.1 Optimization method 

The parameter identification is done by minimizing the discrepancy between a 
reference ECG and the EGG provided by the model with a given set of param- 
eters. More precisely, let us denote by n e N* the number of parameters and 
by 6* e M" the vector of parameters we are looking for in a subset I of M" . The 
subset I is given by Ii x • • • x I„ where Ij is an interval where the value 9j is 
assumed to be. The following cost function is minimized 

J{e) = [ \Vi{t) - yi,ret(t)P + \Vuit) - VuMtW + \Vinit) - Vm.ram'' dt 
Jo 

with respect to 6 £ I where Vi, Vu and Vm are the three Einthoven leads 
given by the simulation for the value 9 of the parameters and V[_rcf , Vu^^ef and 
Viii.rcf are the Einthoven leads of the reference EGG. This functional will be 
approximated by 

N 

J{9) =StY, [\VI,^ - Vl,.ef,,|' + \Vn,^ ~ Vu.,ofaf + \Vul^^ - Vui,rcf ,^\''] 
i=\ 

where Wi is the numerical approximation of W{ti) for W = V\, V\\, V\\\, V[,ref, 
V[i,rcf and Viii^rcf- This cost function only uses information coming from the 
three Einthoven leads, numerical tests have been also run with the twelve stan- 
dard leads and give similar results (at least for the present test case). 

This optimization problem is solved by an evolutionary algorithm (see |lll 
[S]). We briefly indicate its main steps, for the sake of completeness. First, an 
initial population of elements, called individuals, is randomly created. Then, 
the algorithm modifies the population in order to promote the best individuals 
according to the cost function. Gonsidering elements {9i, . . . ,6^^) G I^p 
corresponding to different values of the parameters to identify, the population is 
regenerated Ng times, where Ng corresponds to the number of generations. At 
each generation, J is evaluated for each individual and the population evolves 
from a generation to another following three stochastic principles inspired from 
the darwinian theory of evolution of species 

• selection: promote the individuals whose value by the cost function is 
small, the members which cost function is smaller are preserved in the 
next generation, while those with a high cost function are killed; 

• crossover: create from two individuals two new individuals by doing a 
barycentric combination of them with random and independent coeffi- 
cients, a linear combination of the two selected vectors is done in order to 
create a new member 9 e M" ; 

• mutation: consist of replacing an individual by a new one randomly chosen 
in its neighborhood i.e. a member whose values are closer to the referred 
one in the sense of an Euclidean norm. In our case, the amplitude of the 
mutation get smaller when the number of generations increases. 

At each generation, a one-elitism principle is added in order to make sure to 
keep the best individual of the previous population. 

Among its advantages, the genetic algorithm can easily be run in parallel. 
Moreover, unlike deterministic descent methods which require the computation 
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of the gradient of the cost function, it can be very easily implemented. Its 
main flaw is to require a large number of evaluations of the direct problem, 
and since the initial population is chosen randomly, the minimization process 
has to be run several times. To reduce the number of exact evaluations, we 
have considered the approximate genetic algorithm based on a surrogate model 
strategy (see [Sj). The idea is to approximate the cost function for a part of 
the population, taking advantage of the growing database of exact evaluations. 
The approximation is done by interpolation with Radial Basis Functions. The 
maximum number N^^ of total exact evaluations is flxed and the number of exact 
evaluations decreases at each generation. For instance, in the case of an initial 
population of Np — 80 members and a flxed number of Ng = 15 generations, we 
impose the maximum total number of exact evaluations N^x = 600; typically we 
impose all evaluations to be exact during the flrst 4 or 5 generations, then the 
number of exact evaluations decreases constantly at each generation (see [5 ) . 

Even with the approximate genetic algorithm, solving the optimization prob- 
lem is still very time-consuming. Our strategy is thus to speed up the evaluation 
of the cost function by using the reduced-order model. 

4.2 Identification of four parameters 

In this section, we test our identiflcation strategy for the four parameters Tin, C'mj 
and t^JJ^^,. The reference EGG used in the cost function is the one computed 
in |3] obtained with the full-order model for the values (rin, Cm, An, 7"5^c) = 
(0.8,10-3,200,120). The values of (Ti„, Cm, ^m, r^^oso) are searched for in the 
set 

[0.5,1.5] X [5 X 10"^2 X 10"'''] x [100,300] x [50,150]. 

The key point is to perform the "exact" evaluations required by the optimiza- 
tion algorithm with the reduced-order model. Results presented in Table |2] are 
obtained with 80 POD modes and with the following parameters for the genetic 
algorithm: Np = 120, Ng — 12 and N^x = 850. Just to give an idea, the 
computational time was about 2 days on a PG with 16Go of RAM and using 
six processors Intel Xeon 3.2 GHz. Each one of the 800 time-steps of the "ex- 
act" evaluation requires about 3.5 seconds using the full-order model, while it is 
reduced to about 0.5 second using the reduced model. If the 850 "exact" eval- 
uations were performed with the full-order model, the computational time will 
be unacceptable (more than two weeks). The time needed to compute the POD 
basis is negligible with respect to the time needed by the overall simulation. 

As explained in section |3] the accuracy of the reduced-order model is sat- 
isfactory when the four parameters of interest vary within a reasonable range. 
It is therefore possible to work with a POD basis generated for a single set of 
parameter Oq. This is the approach referred to as Ml in Tables [2] and [3] Never- 
theless, to improve the accuracy, it is possible to use many POD bases computed 
"off-line" for different values of 9 taken in a flnite subset A of /. Next, for a 
given value 9 £ I, the POD basis used for the resolution corresponds to the 
closest value of 9 E A. This is the approach M2 in Tables [2] and |3] 

The relative error is the mean relative error given by the following formula: 

1 / |rm~0.8| [Cm -10-^1 l^m - 200| |rgL-120| \ 
4 V 0.8 10-3 200 120 J ' 
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('''in, Cm, T^^^) 


Relative error (in %) 


Value of J 


Ml with 6*0 = 0^ 


(0.95,9.3 X 10-^ 185,126) 


9.6 


3.05 


Ml with 6*0 = eg 


(0.93,1.05 X 10-3,162,128) 


11.7 


7.4 


M2 with A given by (111 


(0.86, 10-^ 179, 123.5) 


5.2 


2.15 



Table 2: Identification of (Tin, Cm, An, T^fJ^^) (Reference value 
(0.8,10-3,200,120)). 



Here, 6^ = (1, IQ-^, 200, 100), 9^ = (1.3, 1.4 x lO'^, 170, 90) and in the last line 
of the table, we have considered the POD method M2 with A given by 

A= {0.5; 1; 1.5} X {5 X 10--*; 10-3; 1.5 X 10-3; 2 X 10-3} 

X {100; 200; 300} x {50; 100; 150}. ^ ' 

The genetic algorithm is run several times and the results shown corresponds 
to a mean value. On this example, we see that method M2 allows approximately 
to divide the error by 2. The price of this improvement was to compute 108 POD 
bases. Although this computation was done off-line, this approach therefore 
requires a significant computational effort. A larger population has also been 
considered in Table[3]with Np = 300, Ng = 20 and N^x = 1700. Here agam, we 
notice that method M2 allows to improve the identification results. 





A _RV \ 

\iin, ^m, ^m, 'closed 


Relative error (in %) 


Value of J 


Ml with 6*0 = 0^ 


(0.83, 1.02 X 10-3, 184.2, 123.1) 


4.1 


2.2 


Ml with 6*0 = O'f^ 


(0.91,1.09 X 10-3,153,126.4) 


12.3 


7.55 


M2 with A given by ( 


11 


1 


(0.83,1.01 X 10-3,189,123.2) 


3 


2.05 



Table 3: Identification of (Tin, Cm, An, t^^^) (Reference value 
(0.8, 10-3, 200, 120)) with a larger population. 



4.3 Identification of an infarcted zone 



We are interested in estimating the location of an infarcted area, modeled as 
explained in section |3.3[ We proceed as before by minimizing the discrepancy 
between a reference and a simulated EGG, using the genetic algorithm. 

The reference EGG is obtained solving the complete model for an infarction 
area centered at point P as described in section |3.3[ The genetic algorithm is 



run with Np = 80, Ng 
identified are 9 = 
ventricle domain. 



15 and A'e 



600, and the parameters we try to 



{xp.yp, zp). The point P is only searched for in the left 



The results are reported in Figures 12 and 13 The reference EGG (blue line 



in Figure 12 1 is well approximated by the one obtained from the resolution of 



the genetic algorithm (blue line) . The identified infarcted region is actually very 
close to the reference one (Figure 13 1. The solution of the genetic algorithm can 
be improved by including more off-line experiments, as indicated at the end of 
section 13.31 
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Figure 12: Simulated ECG for an infarction located in point P: green 
line represents the simulated reference ECG and blue line the ECG 
corresponding to the infarcted center found with the resolution of a 
genetic algorithm. Black line gives the healthy reference case. 

Time = 80.0 




(a) t = 80ms 



Time = 300.0 




(b) t = 300ms 



Figure 13: Left: transmembrane potential calculated in the reference case 
solved with the full model. Right: transmembrane potential of the solution 
found with the genetic algorithm obtained with the POD. 
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For a different approach to tackle this problem, we refer to |22| . An inter- 
esting possibility to investigate could be to enrich the cost function giving more 
weight to the ST deviation, as suggested in |13j . 

5 Conclusion 

We have presented some preliminary results obtained with a reduced-order 
model of electrophysiology based on the Proper Orthogonal Decomposition 
(POD) method. A well-known difficulty of reduced-order modeling is to identify 
those situations where the reduced-model can be considered as reliable. We do 
not claim that we have answered this difficult question in this paper. Neverthe- 
less, our experiments may suggest some general trends. First, the reduced-order 
model seems quite robust when the ionic current parameters are perturbed and 
when the heart beat increases. This may have interesting applications for long 
time simulations, for example to compute a restitution curve, or to estimate 
the ionic current parameters in an optimization loop. Second, a reduced-order 
model based on a single simulation is totally inadequate to approximate a situa- 
tion with a spatial change in the parameter: this has been shown in the present 
study for an infarcted region; in [4 the same observation was made for the 
initial activation region. A natural strategy consisting of precomputing several 
POD bases with different sets of parameters, or a single POD basis from dif- 
ferent experiments, proved to be satisfactory in the cases we have considered. 
Nevertheless, this solution requires an important off-line effort which makes it 
difficult to apply with more than a few parameters. Alternative strategies have 
therefore to be investigated in future works. 
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